Greatly enhancing the modeling accuracy for distributed parameter systems by 

nonlinear time/space separation 

Hai-Tao Zhang 1 , Chen-Kun Qi 2 , Tao Zhou 3 and Han-Xiong Li 2 

Department of Control Science and Engineering, 
Huazhong University of Science and Technology, Wuhan 430074, PR China 
2 Department of Manufactory Engineering and Engineering Management, 
City University of Hong Kong, Hong Kong SAR, PR China 
3 Department of Modern Physics and Nonlinear Science Center, 
University of Science and Technology of China, Hefei 230026, PR China 
(Dated: February 2, 2008) 

An effective modeling method for nonlinear distributed parameter systems (DPSs) is critical for 
both physical system analysis and industrial engineering. In this paper, we propose a novel DPS 
modeling approach, in which a high-order nonlinear Volterra series is used to separate the time/space 
variables. With almost no additional computational complexity, the modeling accuracy is improved 
more than 20 times in average comparing with the traditional method. 
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Introduction - Most of the physical processes (e.g. ther- 
mal diffusion process P IE IE fl IE IE > thermal radi- 
ation process g|, distributed quantum systems [E ITo| - 
concentration distribution process 

111 El E3, crystal 
growth process PE|, etc.) are nonlinear distributed pa- 
rameter systems (DPSs) with boundary conditions de- 
termined by the system structure. Thus, it is an urgent 
task to design an effective modeling method for nonlin- 
ear DPSs. The key problem in the design of nonlinear- 
DSP modeling method is how to separate the time/space 
variables. Some modeling approaches are previously 
proposed: These include the Karhunen-Loeve (KL ) ap- 
proach P 0, IbH Il5j . the spectrum analysis |l6j. the 
singular value decomposition (SVD) combined with the 
Galerkin's method P[l7|], and so on. Among them, the 
KL approach is the most extensively studied and the 
most widely applied one. In this approach, the output 
T(z, t) is expanded as 

N 

T(z,t) = J2 c i( z )k(t) = C(z)L(t), (1) 

i=l 

where z and t are the space and time variables, re- 
spectively. This operation can be implemented by spa- 
tial basis {ci(z)} combined with time-domain coefficients 
{h(t)}, or time-domain basis {h{t)} combined with spa- 
tial coefficients {ci(z)}. The basis could be Jacobi se- 
ries EH; orthonormal functional series (OFS, such as 
Laguerre series Il8l Il9j . Kautz series [l8|, etc.), spline 
functional series |2J|, trigonometric functional series, or 
some others. However, Nno matter how elaborately the 
basis is designed, the infinite-dimensional nature of DPSs 
does not allow being accurately modeled with small num- 
ber of truncation length N of the basis series. More- 
over, the nonlinear nature of the DPSs will even increase 
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this modeling difficulty. Thus for nonlinear DPSs' mod- 
eling, the extending of N to a sufficiently large number 
is generally required, which would definitely increase the 
computational burden. Consequently, in order to im- 
prove the efficiency of the modeling algorithms, many 
former efforts focused on designing suitable time-domain 
basis {h(t)} or spatial basis \cj{z)\ according to the prior 
knowledge of the DPSs PE|- In addition, some schol- 
ars also presented neural networks to model the non- 
linearities of transitional flows and distributed reacting 
systems based on prop er orthogonal decomposition and 
Galerkin's method |2lL |22| . However, if the prior knowl- 
edge is unavailable or inadequate, the general design 
methods of the bases are very limited so far. On the other 
hand, the conventional finite difference and finite element 
method often lead to very high-order ODEs which are 
inappropriate for dynamical analysis and real-time im- 
plementation [23 . Another conventional approach, spec- 
tral method [l(j, is popularly used for modeling DPSs 
because it may result in very low dimensional ODE sys- 
tems. However, it requires a regular boundary condition 

PH. 

Thus, in this paper, we argue that the linear separa- 
tion is a bottleneck to better modeling performance, and 
to introduce some kinds of nonlinear terms may sharply 
enhance the performance, since they have the capability 
to compensate the residuals of the linear separation. 

The Implement of Nonlinear Space/Time Separation 
- For nonlinear lumping systems, if their dependencies 
on past inputs decrease rapidly enough with time, their 
input /output relationship can be precisely described by 
Volterra series |2E HE HE Ell , which is a gener- 
alization of the convolution description of linear time- 
invariant to time-invariant nonlinear operators. This 
kind of system is called fading memory nonlinear sys- 
tem (FMNS) |2E|, which is well-behaved in the sense 
that it will not exhibit multiple steady-states or other 
related phenomena like chaotic responses. Fortunately, 
most industrial processes are FMNSs. Accordingly, one 
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can naturally extend the concept of Volterra series from 
lumping systems to DPSs by allowing each kernel to con- 
tain both time variable t and space variable z, and then 
design the time/space separation method via the so-call 
distributed Volterra series (see Fig. 1 for the mecha- 
nism of this modeling method). Firstly, the system out- 
put can be represented by 

T(z,t) = h (z) + / hi(z,Ti)u(t — Txjdrx + 



Jo 



h 2 (z, Ti,T 2 )u{t - T\)u{t - T 2 )dT 1 dT 2 H , (2) 



where hi(z, ri,--- , r,) is the ith order distributed 
Volterra kernel. Then denote 4>i (t) as the ith order OFS 
and li(t) = f Q <f)i(t)u(t — r)dr as the ith order OFS filter 
output. Since {<j>i(t)} forms a complete orthonormal set 
in functional space l 2 , each kernel can be approximately 
represented as 



hi(z,n) = y^c fc (z)0fc(ri), 

k=l 

N N 

h 2 (z, Ti,T 2 ) = ^ X! C nm(z)(f>n(Tl)(t) m (T 2 ), (3) 



where Cfc(z) and c nm (z) are spatial coefficients. Then, 
the input /output relationship can be written as (see Eq. 
(1) for comparison) 



T(z, t) = c Q (z) + C(z)L(t) + L T (t)D(z)L(t) + ■ 



(4) 



where L(t) = [h(t) ■ ■ ■ l N (t)] T , C(z) = [a(z) ■ ■ ■ c N (z)], 
and D(z) = [c i:j (z)] NxN . 

To obtain the spatial coefficients, firstly we pre- 
compute all the OFS kernels according to the polynomial 
iterations [3 or the following state equation 



L(t+1) = AL(t)+Bu{t), 



(5) 



where u(t) is the system input, and A and B are pre- 
optimized matrices (see Ref. [2ij for details). Then the 
input / output relationship Eq. (4) is represented by a lin- 
ear regressive form, and these spatial coefficients cq(z), 
C(z), D(z), ■ ■ ■ can be obtained by using the least square 
estimation combined with spline interpolation [2(]]. Fi- 
nally, the model is obtained by synthesizing the OFS ker- 
nels and the spatial coefficients according to Eq. (4). 

Fig. 2 shows the operation details of this modeling 
method. The first order OFS filter is the Laguerre series, 
in which 



Goiq- 1 ) 



-,Gi(r X ) 



q 



i-<r 



(6) 



where a is the time-scaling constant and q~ l is the one 
step backward shifting operator (i.e. g -1 w(£) = u(t — 1)). 
The second order OFS filter is the Kautz Series, in which 
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FIG. 1: The sketch map of OFS- Volterra modeling for non- 
linear DPS. 
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FIG. 2: Operation details of OFS- Volterra modeling. 



Go(q~ 1 ) and Gi(q~ 1 ) are second order transfer functions. 
Analogically, Heuberger et dl. [3(| introduced the higher 
order OFS model. As the order increases, OFS model 
can handle more complex dynamics. 

Numerical Results - Consider a long, thin rod in a re- 
actor as shown in Fig. 3. The reactor is fed with pure 
species A and a zeroth order exothermic catalytic reac- 
tion of the form A — > B takes place in the rod. Since 
the reaction is exothermic, a cooling medium that is in 



cooling me dium 



A-> B 



A 



B 



FIG. 3: The sketch map of catalytic rod. 
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FIG. 4: System output. 
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FIG. 5: Modeling error of the traditional method. 



contact with the rod is used for cooling. Assume the 
density, heat capacity, conductivity and temperature are 
all constant, and species A is superfluous in the furnace, 
then the mathematical model which describes the spa- 
tiotemporal evolution of the rod temperature consists of 
the following parabolic partial differential equation: 



dT _ d 2 T 
dt dz 2 



+ p T e~— -fae-i + (3 u (b(z)u(t) - T) , (7) 



e 2 (z,t) 
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FIG. 7: IAEs of the traditional method (left) and the present 
method (right). 




FIG. 6: Modeling error of the present method. 



FIG. 8: ITAEs of the traditional method (left) and the present 
method (right). 



which subjects to the Dirichlet boundary conditions: 

r(0,t) = 0, T(vr,<) = 0, T(z,0) = 0, (8) 

where T(z,t), b(z), Pt, flu, 7, and u denote the tem- 
perature in the reactor (output), the actuator distribu- 
tion function, the heat of reaction, the heat transfer co- 
efficient, the activation energy, and the temperature of 
the cooling medium (input), respectively. Here we set 
f3 T = 50.0, (3 U = 2.0, and 7 = 4.0. In the numerical cal- 
culation, without loss of generality, we set b{z) = 1, and 
u(t) = [1.4,1.4,1.4,1.4]. The order of Volterra series is 
two, the OFS is chosen as one-order Laguerre series 0] 
with a = 0.6, and the truncation length is given N = 4. 

The system output is shown in Fig. 4. Denote by 
e(z, t) the modeling error, that is, the difference between 
system output and modeling result at the point (z,t). 
Fig. 5 and Fig. 6 exhibit the modeling errors of the tradi- 
tional and the present methods, respectively. Clearly, the 
method proposed here has remarkably smaller error than 
that of the traditional one. To provide a vivid contrast 
between these two methods, we calculate the integral of 
absolute error (IAE, J \e(z,t)\dz) and time-weighted ab- 
solute errors (ITAE, J t\e(z,t)\dt), which are two stan- 
dard error indexes to evaluate modeling performances of 
DPS and can be considered as the modeling errors along 
the temporal dimension t and the spatial dimension z. As 
are shown in Fig. 7 and Fig. 8, both the IAE and ITAE 
of the present method is reduced by > 20 times com- 
paring with those of the traditional one, which strongly 
demonstrates the advantage of the present method. Ac- 
tually, to obtain the shapes of IAE and ITAE, one can 
cut the error surfaces of Figs. 5 and 6 along the tem- 
poral coordinate t and the spatial dimension z. In addi- 



4 



1.14 
1.11 
1.08 
1.05 
1.02 
0.99 





3 4 5 6 7 



3 4 5 6 7 8 9 



FIG. 9: Ratios of the consumed time tz/tj (left) and the av- 
erage of absolute modeling error |ei|/|ea| (right) vs the trun- 
cation length N. The subscripts 1 and 2 denote the cases of 
the traditional and the present methods, respectively. The 
CPU time by using traditional method for N 6 [3, 9] is in the 
interval [100s, 160s]. All the numerical calculations are imple- 
mented by using a personal computer with CPU: 1.8G, RAM: 
256M, OS: Windows XP, and software platform: MATLAB 
6.5. 
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FIG. 10: Modeling methodology for DPSs. 



tion, we calculate the average of absolute modeling error 
Tj^uSi 1 1 \e(z,t)\dzdt. From Fig. 9, it is found that, 
in comparison with the traditional method, the model- 
ing accuracy of the present one is enhanced by 14-32 
times with less than 15% increase of the consumed time. 
It should be note that the modeling accuracy would in- 
crease along with the increase of the Volterra series order 
N v , however, the computational complexity will increase 
too. For lumping systems, this fact has been proved, and 
for DPSs, this fact is also validated by experimental re- 
sults [24| . So a tradeoff between the modeling accuracy 
and the computational complexity must be made. That 
is why here we set the order N v = 2. 

Conclusion and Discussion - Modeling method for 
nonlinear DPS plays an important role in physical sys- 
tem analysis and industrial engineering. Unfortunately, 
there exits two essential difficulties in this issue, a) 
distributed nature due to time-space coupled, which 
causes different temperature responses at different loca- 



tions; b) nonlinear complexity from varying working 
point - different dynamics obtained even at the same lo- 
cation for a large change of working points. Owing to 
these difficulties, previous modeling methods via linear 
time/spatial separation techniques (e.g. KL approach, 
spectrum analysis, SVD-Galerkin technique, etc.) can 
not yield satisfying modeling performance, especially for 
DPSs with severe nonlinearity. The modeling error is 
caused by the nonlinear residue of the linear time/space 
separation. Thus, it is naturally to expect that a non- 
linear time/space separation method may yield better 
modeling performance. 

To validate this supposition, we design a novel model- 
ing method by extending the concept of lumping Volterra 
series to the distributed scenario. As shown in Fig. 10, 
the nonlinear DPS is first decomposed into kernels, upon 
which the time-space separation is carried out. These 
two decompositions will gradually separate the com- 
plexity and provide a better modeling platform. The 
time/space separation will be handled by a novel KL 
Volterra method instead of the conventional KL method, 
the time-domain complexity by the OFS-based learning, 
and the spatial complexity by the curve fitting techniques 
(e.g. spline interpolation) or intelligent learning algo- 
rithms (e.g. neural network, fuzzy system, etc.). This 
novel method is applied on a benchmark nonlinear DPS 
of industrial process, a catalytic rod. It is found that the 
modeling accuracy is improved by more than 20 times 
in average comparing with the traditional method, with 
almost no additional computational complexity. The un- 
derlying reason may be that the high order Volterra ker- 
nel can compensate the residuals of the linear separa- 
tion. In addition, we have applied this method to an- 
other two benchmark nonlinear DPSs, a rapid thermal 
chemical vapor deposition process Q, and a Czochral- 
ski crystal growth process [lj . The corresponding results 
also strongly suggest that the nonlinear time/space sep- 
aration can greatly enhance the modeling accuracy. 

Although its superority has been demonstrated, the 
KL Volterra method is just a first attempt aiming at the 
motivation of nonlinear time/space separation. Thanks 
to its excellent modeling efficiency, this novel method is 
definitely a promising one for both physical system anal- 
ysis and industrial engineering. We believe that our work 
can enlighten the readers on this interesting subject. 
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